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Abstract 


Differential equations of fractional order have been the focus of many stud- 
ies due to their frequent appearance in various applications in fluid me- 
chanics, biology, physics, and engineering. In general, it is not easy to 
derive the analytical solutions to most of these equations. Therefore, it is 
vital to develop some reliable and efficient techniques to solve fractional dif- 
ferential equations. A numerical method for solving fractional differential 
equations is proposed in this paper. The method is based on a hybrid of 
Block-pulse and orthonormal Bernstein functions. Convergence analysis is 
given, and numerical examples are introduced to illustrate the effectiveness 
and simplicity of the method. 
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1 Introduction 


Fractional calculus is one of the most popular calculus types having a 
vast range of applications in many different scientific and engineering dis- 
ciplines. The order of the derivatives in the fractional calculus might 
be any real number that separates the fractional calculus from the ordi- 
nary calculus where the order of derivatives are allowed to be only nat- 
ural numbers. Fractional calculus is a highly efficient and useful tool in 
the modeling of many sorts of scientific phenomena including image pro- 
cessing, earthquake engineering, biomedical engineering, computational fluid 
mechanics, and physics. Fundamental concepts of fractional calculus and 
its applications to different research areas can be seen in the references 
[3, 4, 2, 7, 6, 5, 9, 12, 10, 16, 19, 20, 24, 29, 31, 32, 34], amongst many others. 
In recent years, many different basis functions have been used for solving 
fractional equations, such as operational matrix of Chebyshev polynomials, 
Block-pulse functions (BPFs), hybrid Bernoulli and Block-pulse functions, 
and hybrid Legendre [1, 8, 15, 23, 27, 26, 30, 33]. 

In this paper, we employ hybrid functions consisting of a combination of 
BPFs with orthonormal Bernstein polynomials (OBPs) to find the numerical 
solution of the fractional equation 


F(a,y(x), D“ y(a),...,D° y(a)) =0, « € [0,X], a > 0, t=1,...,k, 
with boundary or supplementary conditions 


Zily(Ei).y (E),---,y (Ei) =di, 4=0,1,...,p, (2) 


where 0 < p < max{aq;,i = 1,...,k} <p+1, & © [0,X], i =0,...,p and 
Z;,i = 0,...,p, are independent linear combinations of y(€;), y (&),-..,y) (&) 
and y € L?/0, X]. It should be noted that F can be nonlinear in general. 
The fractional equations are used in a variety of fields including con- 
tinuum mechanics, potential theory, geophysics, electricity and magnetism, 
antenna synthesis, communication theory, mathematical economics, popu- 
lation genetics, the particle transport problem in astrophysics, and reactor 
theory. For the most part, it can be challenging to deduce the analytical 
solutions to most of these equations. As a result, it is crucial to develop 
some reliable and efficient ways to solve fractional differential equations. In 
this paper, we introduce a new method based on a hybrid of Block-pulse 
and orthonormal Bernstein functions. The present study’s primary goal is 
to present and analyze a new stable algorithm for the numerical solution of 
the fractional differential equation based upon the hybrid of Block-pulse and 
orthonormal Bernstein functions (HBB method). To solve this problem, we 
propose a mathematical formulation that takes advantage of the orthogo- 
nality of BPFs and the OBPs properties to reduce the fractional differential 
equation to an algebraic system. Convergence analysis of the method is 
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discussed theoretically and numerically. Moreover, the applicability of the 
method is examined by solving some fractional differential equations and the 
Basset equation, which describes the unsteady motion of a sphere immersed 
in a Stokes fluid. 

The outline of this paper is as follows: In section 2, we briefly overview 
some basic definitions for fractional calculus. In section 3, hybrid functions 
are introduced; therefore we approximate functions by using hybrid func- 
tions. In section 4, we present useful properties of hybrid functions such 
as coefficient matrix and operational matrix of derivative to solve fractional 
differential equations. In section 7, we present estimates for the error of the 
best approximation of smooth functions by hybrid functions. In section 6, 
we apply the proposed method to find an approximate solution to fractional 
differential equations. Finally, numerical examples are presented to show the 
effectiveness of the proposed method in section 7. 


2 Fractional calculus 


We give some basic definitions and properties of the fractional calculus theory, 
which are used further in this paper. 


Definition 1 (see [32]). The Riemann-—Liouville fractional integral operator 
of order a > 0 is defined as 


Fee a fe -) 70d = fe), oS 0,450) 
J° f(x) = f(z). (4) 
It has the following property: 


Whar cal _ Cy + by ito 


. sa 
I(y+1+a) : 


Definition 2 (see [32]). The Caputo definition of fractional derivative oper- 
ator is given by 
1 


D*f(e) =J"-°D" fe) = 


[ @-9r oa, 
) Jo 
where m—1<a<m, mEN, x > 0. For the Caputo derivative, we have 


D°C=0 (C is a constant), 


. 5 
TG+1) j-a for 7 € NUO and j > [a] or f EN. (5) 


Daf" for 7 € NU0 and j < [a], 
T(j+l—a) 
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3 Hybrid of Block-pulse and orthonormal Bernstein 
functions 


The orthonormal set of hybrids h;;(x), 7 = 1,2,...,m andi = 0,1,...,n, 
where j is the order for BPFs, 7 is the order for OBPs, and z is the normalized 
time, is defined on the interval [0, X) as 


(6) 


hte Vm bin(ma —y) if SEY <g < XS, 
Mee 
0 otherwise, 


where y; = X(1—j) and bj,,i =0,...,n are OBPs defined on [0, X]. Then 
H (a) = [hin (x), sey hin(2), sey hjo(2), sey Rige ha), sey Amo(x), sey hienlayy” 


=|Fi (a), 2c Aye): ae) 


is the m(n + 1) hybrid function vector. Here, the Bernstein polynomials 
(BPs) of the nth degree are defined on the interval [0, X] as 


1 
Bin(x) = ro (") “(xX —a2)"", i=0,1,2,...,n, 


by using expansion of (X — «)"~', we have 


where 


According to [25], BPs form a complete basis over the interval [0, X]. 


The explicit representation of the OBPs, denoted by bj,(x) here, was 
recently discovered by analyzing the resulting orthonormal polynomials after 
applying the Gram—Schmidt process; see [13] 


a 


bina) = (VB DAD EX—ay" See (PE 8) (ars 


k=0 


In addition, it can be written in terms of the non-orthonormal Bernstein basis 
functions as 
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bin(x) = (\/2(n — i) 4 DS 1)" en ’ (:) 


T van oon i—-kn—k(X), 
7 i-k 
for i =0,1,...,n. By using the Taylor expansion, };,,(a) can be represented 
as [21] 
bin = ZiT, (2), 7=0,...,n, 
where T;,(x) = [1,2,27,...,2"]", Z; is a row vector of Taylor coefficients as 


min{i,j} 
(Zi)3=V2n—-)+1 So pij—nBix, J =0,...5n, 


k=max{0,j—n+i} 


where 


,_,f2nt1l—-it+g i ; ; 
Bij = (-1) ( j ee j=0,...,%. 


Thus (2) = [bon(z), bin(x),---,;bnn(x)|* can be defined by 
b(x) = ZT, (2), (7) 
where Z is an (n+ 1) x (n+ 1) matrix whose ith row is Z;,1 =0,...,n. 
Now, from (7) and (4), we have 
H;(x) = /mZT, (mx — 4;), (8) 
where 
Tp(ma — 74) = [1,ma — 74, (ma —4)®, +, (ma =)", 
by using the binomial expansion of (ma — 7;)", (8) can be expressed as 


H;(x) = ZA;Tn (2), 


where 
1 0 0 0 tee 0 0 
V4 m 0 0 tee 0 0 
s=vm} Oa)? Gam (2)? 0 a 0 0 


(6) 9)” (1) ag) (5) Cg)" Pm? (5) (ag)? Ga) (a)m™— (7) me” 
Therefore the hybrid function vector H(a) can be defined by 
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A(x) =VT(«), (9) 
where T(x) = [Tn(x),Tn(x),---,Tn(a)]? is an m(n +1) vector, and T = 


diag[T1,T2,..., Pm] is an m(n+1) x m(n + 1) coefficient matrix of the (n + 
1) x (n+ 1) coefficient submatrix I; = ZA;. 


3.1 Function approximation 


Suppose that Q = L?[0, X] and 
Lad yeep hing 4 yO ang Rares Dinas ee nin} cQ 
are the set of hybrid functions and that 


Qmn = spant{hio,...,hin,-.-,hjo,---,hyns--+,Rmo,---; Amn}. 


Let fmn be the best approximation of an arbitrary function f € L?(Q) 
out of Qmn, that is, 


IIf— fmnll2 Sil fg lle for all g € Qmn. 


Moreover, since f,,,, € Qmn, there exist unique coefficients C19, C11,---,€mn 
such that 


fle) = ful) = Yo eyshy(@) =D) CH (2) = C7H(2) = H™(@)C, 


where C and H(a) are m(n + 1) vectors as 


fi T 
C= [C10s=' +29 Clvigs +5 CfOjno 4 Cpa ses Cmte 5 Gan = [C1,...,Cj,...,Cm] ’ 


HS ilgyies Nines re lsly oes Ritas veg enbineca nel Aiea aval’: 


and the hybrid function coefficients cj; are obtained by 


Since the set of hybrid functions is a complete orthonormal system in 2, then 


x 
|Ajill2 = ip hole )h,(ejdz = 1. 


Therefore 
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x 
oj = (Fs hyi) = | P(a)hyi(w)de. (10) 


4 Operational matrix of differentiation 


We construct the operational matrices of differentiation for H(x). First, we 
consider integer order derivatives of H(x). Hence 
dH (x) 
dx 


=DOH(2), 


where D™) is the m(n + 1) x m(n + 1) operational matrix of first derivative 
of Hybrid function. From (9), we have 


dH(x) _ aT (a) 


ae a =TQT(x) =TQr-!A(a), 


where Q = diag[Q,Q,...,Q] is an m(n+1) x m(n +1) matrix that consists 
the m submatrix Q(n41)x(n+41) a8 


000---00 
100---00 
g=| 020-00 


000---n0 


(n+1) x (n+1) 
Therefore we have Ps 
pb” =Trer-!. (11) 
In general, we obtain 
d"H 
© _ (pMyrH(@), 
dx” 


where n € N and the superscript, in D), denotes the matrix powers. Thus 


D®™ =(D™)", n=1,2 


i 


For the Caputo fractional derivative introduced in section 2, we denote frac- 
tional order derivative of H(a) by 


where D(® is an m(n +1) x m(n +1) operational matrix of fractional order 
derivative of Hybrid function. From (9), we have 
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d° H(x) pero 


dx@ dx 


>) 


where 


d°T (x) = Se EY Ee) aC 
dx —* dx® * daw ?"""’ 


and by using (3), we can specify the elements of Fr) as 


da _ J 9, 0<k< fal, 
dz® ree, fa] <k<n. 
We can also approximate «*~° for k = [a],...,n by hybrid functions as 


gk = S > prgshge(2) = PP H(z), 


where P; is an m(n + 1) vector that we can obtain by (10) and (9) as 


x Xi 
Pei = iE a*-°h,;,(a)da =f; cee gh-atign 
Therefore 
Do = TWP, (12) 


where © = diag[V,¥,...,V] is an m(n +1) x m(n +1) matrix with m 
submatrix W(n41)x(m41) as 


T(fa] + 1) T(n+1) 


U = diag(0,... oe 
AAO Dien One cea ay Tek ba) 


F 


and P = [P,P,..., P|" is the m(n +1) x m(n +1) matrix with the (n+1) x 
m(n +1) submatrix P, defined as 


Pfa]10 Pfaj11 Pfaji2 °** Pia]m(n—-1) Pfa|mn 
P= : A x ‘i 


Pnio Pnil Pni2 *** Pnm(n-1) Pnmn 


5 Convergence analysis 


The purpose of this section is to obtain an estimate of the error norm of 
the best approximation of a smooth function of two variables, on a certain 
domain [0,X], by a bivariate polynomial. This estimate will be used for 
comparisons in section 6, when analyzing the error of the numerical results. 
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We assume that f = pet f; is a sufficiently smooth function on [0,X] 
and that p = yA p; is the interpolating polynomial to f at points 
xj,t = 0,1,...,n, which are the roots of the (n+1)-degree shifted Chebyshev 
polynomial on [0, X]. Then we have [17] 


7 m 7 m rrr) n 
f(z) — P(e) = i - Pi) = IIe x)), 


! 
j=l j=l (n+)! i= 


where € € [Ij = [ews ix), Therefore 


. (n+1) ee: _ @e) 
— < : a 
| f(x) — P(x) Is y a eee os Ve uo) 
We assume that there is a real number y such that 
max | f,"* (2) |< 4. (14) 


By replacing (14) in (13) and taking into account the estimates for Chebyshev 
interpolation nodes [28], we obtain 


Xr 
| f(x) _ P(z) I< TpantImn(n + Yr 


(15) 


With the help of (15) we obtain the following result. 


Theorem 1. Let fimn(x) = C7 H(z) be the hybrid functions expansion of 
the sufficiently smooth real function f in Q, where 


C = [eip(z),.++5Cinl2),...,070(@), 2. tpn (@)y< ++ 5 Cn0(Z), <<05 Creel @)I" 5 
and 
x 
c= fo Ble\hp (eae. 
0 
Then, there exists a real number 7’ such that 


xnrri 
= <o/ . 
\|f fmnll2 ae | Q2nt+1mn(n + 1)! 


(16) 


Proof. Let Qmn be the space of bivariate polynomials of degree < n. As 
shown in section 3, fmn is the best approximation of f in Q»,, that is, 


where g is any arbitrary polynomial in 2,,,,. In particular, we have 
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x 
Ilf — fmnll2dx < 7 | A) — Ppa (at) de, (17) 


where p,,,, is the interpolating polynomial of f. From (15) and (17), we 
obtain 


J x Xn. ‘ xynrti i 
= < dx = X . 


From (18), we conclude that (16) is valid with 


Theorem 2. Suppose f € L?[0,1],n —1<a<n. Then 


” e 1 xrtri 
|| D°f — D* finn \|S : 


; 1 
= Tin—atl)” 22r+1mn(n +1)! ay) 


Proof. By using (18) and [11], we have 


|| D°f — D* fn ||? =|] P°~°(D" f — D” finn) II? 


1 mn nm 
=|| xrito—nT'(n — a) (D f D tram) ||? 
1 2 n n 2 
< ‘ma Sar@ow). |) a a Fra ll 
1 
< (hho = Dy” Il f —fmn |l? . 


From (18), we get (19). 


6 Solving fractional differential equations 


In this section, in order to show the high importance of the proposed method, 
we apply it to solve fractional differential equation (1) and (2). In order to 
use hybrid functions for this problem, we approximate by hybrid functions 


as 
mon 


y(a) ~ S- S > cjihji(a) = CTH(), (20) 
j=l i=0 
where C = [ci0,---;€iny+++3€m0;+++;€mn|~ is an unknown vector. Using (11) 


and (12), we have 


D%y(x) ~ CT D% H(z) ~ CTD) H(z), j =1,...,k. 
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By substituting these equations in (1) and (2), we get 
F(a2,C? H(x), C?D() H(a),...,C?D(*) H(«)) = 0, (21) 


To find the solution y, we first collocate (21) at m(n + 1) — (p +1) points. 
For suitable collocation points, we use 


21-1 


= Sai eats bem +1) - @+). 


Xi 
These equations together with (22) generate m(n + 1) algebraic equations, 
which can be solved to find c;;,j = 1,...,m,i = 0,...,n. Consequently, 
the approximate solution of the unknown function y, given in (20), can be 
calculated. 


7 Numerical examples 


In this section, numerical results of some examples are presented. The abso- 
lute errors of this method are compared with those of the existing methods 
reported in [35, 14, 22]. The computations associated with these examples 
were performed using Mathematica 9.0. 


Example 1. Consider the fractional equation 


ViD?y(x) + e*y(a) = m2 +e"), O<a<l, (23) 


Ti 


with the following initial conditions y(0) = 0. The exact solution of this 
problem is y(x) = x. The numerical method developed in section 6 is applied 
to (23) for m = n = 2, we approximate solution as 


y(&) & ciohio + c11hi1 + Crahi2 + c2ooh29 + catha1 + cozhe2 = CT H(x). 


Here, we have 


3.16228 —12.6491 12.6491 0 0 0 
—2.44949 29.3939 —48.9898 0 0 0 
r- 1.41421 —22.6274 56.5685 0 0 0 
0 0 0 12.6491 —25.2982 12.6491 : 
0 0 0 —29.3939 78.3837 —48.9898 
0 0 0 26.8701 —79.196 56.5685 
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0 0 0 0 0 0 
0 1.12838 0 0 0 0 
yp — 0 0 1.50451 0 0 0 
~ 10 0 0 0 0 0 ? 
0 0 0 0 1.12838 0 
0 0 0 0 0 1.50451 
2.38514 0. 0.666667 0.672262 0.430222 0.236893 
0.170367 0.263932 0.161905 0.415476 0.383286 0.235524 
p-— 0.0283945 0.0879772 0.084127 0.262718 0.331263 0.235838 
2.38514 0. 0.666667 0.672262 0.430222 0.236893 | ’ 


0.170367 0.263932 0.161905 0.415476 0.383286 0.235524 
0.0283945 0.0879772 0.084127 0.262718 0.331263 0.235838 


1.89128 —2.09283 —0.709874 —0.930387 0.8335 1.12652 

3.55781 2.26954 —0.830648 —5.58347 —11.7032 —9.57079 

D:- —1.93327 0.748753 3.02605 11.7513 18.4068 14.0581 
—4,32293 —5.85992 —3.02074 —6.86047 —4.63714 —2.23512 

12.9755 16.8594 8.1193 17.3836 9.48449 3.44875 

—12.8079 —16.0982 —7.30845 —14.7689 —6.05866 —0.975538 


BR 


Finally, the obtained results are reported in Table 1. It shows that our 


method has better accuracy with the piecewise constant orthogonal function 
developed in [35]. 


Example 2. Consider the following fractional differential equation: 


Di y(a) + x3 y(a) = ron” +23, 0<2<1, 
with the following initial conditions 
y(0) = 0. 
We know that the exact solution is y(z) = x. The maximum errors for 


different values of x are listed in Table 1. It is shown that, we obtain the 
error of order 107! for m = 2,n = 2. In the Haar wavelets method [14], the 
error is 0.0014 obtained for m = 64, where m is the order of Haar wavelets. 
Clearly, the results obtained by our method are better than those in [14] in 
terms of accuracy if the exact solution is sufficiently smooth. Therefore, the 
HBB method is a valid method in solving fractional differential equations. 


Table 1: Numerical results for Examples 1 and 2 for m = 2 and n = 2. 
Aa 0.1 0.3 0.5 0.7 0.9 
Example 1 3.46045x 10-17 5.55112 10-7” 3.21965x 10" 233147x 10°" 1.77636x 10” 
Example 2. 1.38778 x 10-17) 5.55112 x 10-17) 1.5099 x 107" 1.55431 x 10-15 1.27676 x 10-' 


Example 3. Consider the fractional equation 
2-a 


D“y(«) _ —y(a) + x + T3—a) 
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with the following initial conditions y(0) = 0. The accurate solution, in 
this case, is given by y(2) = 2”. For a = 0.25,0.75 and m = n = 2, the 
approximation solution for (24) is obtained by using the HBB method. The 
absolute error of HBB method and fast wavelet collocation method (FWCM) 
reported in [22] is shown in Table 2 when a = 0.25, a = 0.75, and m =n = 2. 


The logarithm of absolute error for different values of m and n is shown in 
Figures 1 and 2. Table 2 and Figures 1 and 2 confirm that the HBB method 
approximates the solution of fractional differential equation uniformly. 


n=2,m=1 
4.x 10-181 . ) ) ) ‘| 
5.x 10-16 


1.x 10°76 
5.x10°174 


Log,)(Error) 


1.x 10-17 


107135 4 
107144 4 
10715 5 4 


10776 


Log,,(Error) 


107-17 L 


10-78 L 


10°19 L 


Figure 1: The curves of the logarithm of the absolute errors of Example 3 for a = 0.25 
and different value of m and n. 
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n=3,m=1 


10-14 L 


40-15 L 


Log,,(Error) 


40-16 I 


Log, (Error) 
3 
i 
= 


407-76 L 4 


Figure 2: The curves of the logarithm of the absolute errors of Example 3 for a = 0.75 
and different value of m and n. 


Table 2: The errors of Example 3 compared to [22 
x a = 0.25 a= 0.75 
Current Method —FWC method in [22] Current Method FWC method in [22] 

0.03125 3.46945 x 10-18 (0.05958 x 10-2 3.95517 x 10- © (0.33012 x 10-8 
0.09375 2.77556 x 10-17 ~— 0.04301 x 107” 4.02456 x 10-'® 0.34730 x 10718 
0.18750 8.32667 x 10-17 0.05390 x 10-12 3.88578 x 10-!6 0.02296 x 10-38 
0.28125 1.38778 x 10-!7 0.02446 x 10-!2 3.88578 x 10-16 0.19817 x 10-38 
0.37500 2.22045 x 10716 0.00394 x 107” 4.16334 x 107'® 0.29393 x 10718 
0.46875 6.15064 x 10-!4 0.06483 x 10-12 5.55112 x 10716 0.22148 x 1078 
0.56250 5.77316 x 10-4 0.05401 x 10-12 6.66134 x 107 0.19428 x 10-38 
0.65625 5.54001 x 107-'4 0.00982 x 107!” 8.88178 x 107! 0.17097 x 10718 
0.75000 5.42899 x 10-4 0.01798 x 10-12 1.33227 x 10-15 0.09103 x 10-18 
0.84375 5.45127 x 107-'4 0.03186 x 107” 1.77636 x 107! 0.11990 x 1078 
0.93750 6.25614 x 10-'4 0.20128 x 10-2 2.47315 x 10-' .27533 x 107-8 
CPU time 3.65s 2.578 
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Example 4 (see [18]). Consider the following fractional Relaxation problem: 


Diy(x) =—y(z), O<2<l, (25) 
with initial condition y(0) = 1. 
The exact solution of this equation is y(~) = E,(—a#°). Here E,(2) is 


called the Mittag—Lefler function, 


Table 3 shows a comparison of our method and the Taylor matrix method 
presented in [18]. Clearly, the results obtained by the HBB method are in 
agreement with other mentioned numerical method and in total this approach 
has high accuracy. 


Table 3: Comparison of our method for m = 3 and n = 8, and the Taylor matrix 
method in [18] of Example 4 


x Current Method  Ref.[18] 


0.0 3.4685 x 10-1 ~——-0.0000 

0.1 1.7745 x 10718 0.1000 x 10~° 
0.2 6.3421x 10729 0.9000 x 1079 
0.3 5.3461 x 10-9 0.4000 x 1079 
0.4 6.4326x 10729 0.2000 x 1079 
0.5 7.4902 x 10-9 0.1000 x 10-8 
0.6 4.3721 x 1078 0.6000 x 10-8 
0.7 3.5684 x 1077 0.2400 x 1077 
0.8 9.2341 x 1077 0.8400 x 1077 
0.9 8.4563 x 10-7 0.2480 x 10-8 
1.0 5.3475 x 107-7 0.6740 x 10-8 


8 Conclusion 


In this paper, we obtained operational matrices of the fractional derivatives 
of a combination of OBPs and BPFs. Then by using these matrices, we 
reduced the multi-order fractional differential equations to a system of alge- 
braic equations that can be solved easily. The convergence analysis of the 
HBB bases is given in section 7. The numerical solutions obtained using 
the suggested method showed that numerical solutions are in very good co- 
incidence with the exact solution. For future research, we will apply this 
numerical method for solving nonlinear fractional integro-differential equa- 
tions, nonlinear Volterra—Fredholm—Hammerstein integral equations, and so 
on. 
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